Multitemporal Polarimetic SAR Change Detection

Introduction

The Docker container mort/sardocker provides easy access to Python scripts for analysis of polarimetric synthetic aperture radar (polSAR) imagery. These scripts are described in Canty(2014), Image analysis, Classification and Change Detection in Remote Sensing, 3rd Revised Ed. In addition to scripts for polSAR change detection, the container encapsulates the command line interface of the ASF MapReady software for terrain correction and geocoding of SAR images. The user interacts with the software in an IPython notebook served from within the Docker container.

Software Installation

On 64-bit Ubuntu Linux:

  1. install Docker
  2. In a terminal, run the command
    sudo docker run -d -p 433:8888 --name=sar -v your_image_directory:/sar/imagery mort/sardocker
    where your_image_directory is the path to your SAR data.
  3. Point your browser to
    localhost:433
  4. Click on this tutorial, or open a new notebook with New/Python 2.

On WIndows (or Mac), install boot2docker, share your image directory with VirtualBox and proceed from step 2. above.

Here is a listing of the main directory /sar in the container. It contains the various Python and bash scripts required for preprocessing and change detection:

In [12]:
!ls -l /sar
total 6548
-rw-rw-r--  1 root root   11115 Oct  4 15:26 dispms.py
-rw-rw-r--  1 root root    7749 Oct 15 08:32 enlml.py
drwx------ 27 1000 1000    4096 Oct 12 10:25 imagery
-rw-rw-r--  1 root root    4408 Oct 15 08:00 ingest.py
-rwxr-xr-x  1 root root    8072 Oct 12 19:37 libprov_means.so
-rwxrw-r--  1 root root    2923 Oct 15 09:13 mapready.sh
-rw-rw-r--  1 root root    8591 Oct 15 09:42 omnibus.py
-rwxrw-r--  1 root root    1139 Oct 15 09:29 omnibus.sh
-rw-rw-r--  1 root root     693 May 17 14:58 prov_means.c
-rw-rw-r--  1 root root     841 Oct 10 13:39 radarsat2quadpol.template
-rw-rw-r--  1 root root    8636 Jun 10 12:55 register.py
-rw-rw-r--  1 root root     828 Oct 15 08:32 terrasarxdualpol.template
-rw-rw-r--  1 root root 3891749 Oct 15 12:55 tutorial.ipynb
-rw-rw-r--  1 root root      50 Jun 10 11:14 utm.prj
-rw-rw-r--  1 root root 2699882 Sep 22 09:12 vortrag.ipynb
-rw-rw-r--  1 root root   10262 Oct  6 12:34 wishart.py
-rwxrw-r--  1 root root     812 Oct  5 07:20 wishart.sh

The /sar/imagery directory contains the polarimetric SAR data and is shared with the host. In the present case there are 12 Radarsat-2 quadpol images in SLC (single-look complex) format along with a dem (digital elevation model). Acquistion times range from May 25, 2009 (20090525) to October 11, 2010 (20101011):

In [11]:
!ls -l /sar/imagery | grep "_SLC$"
drwx------ 4 1000 1000 4096 Oct 13 11:54 RS2_OK5491_PK71074_DK68879_FQ21_20090525_172447_HH_VV_HV_VH_SLC
drwx------ 4 1000 1000 4096 Oct 11 15:25 RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC
drwx------ 4 1000 1000 4096 Oct 11 13:16 RS2_OK5491_PK71074_DK68879_FQ21_20100520_172458_HH_VV_HV_VH_SLC
drwx------ 4 1000 1000 4096 Oct 11 14:11 RS2_OK5491_PK71074_DK68879_FQ21_20100707_172459_HH_VV_HV_VH_SLC
drwx------ 4 1000 1000 4096 Oct 11 14:12 RS2_OK5491_PK71074_DK68879_FQ21_20100731_172501_HH_VV_HV_VH_SLC
drwx------ 4 1000 1000 4096 Oct 11 14:13 RS2_OK5491_PK71074_DK68879_FQ21_20100824_172503_HH_VV_HV_VH_SLC
drwx------ 4 1000 1000 4096 Oct 11 14:14 RS2_OK5491_PK71074_DK68879_FQ21_20101011_172507_HH_VV_HV_VH_SLC
drwx------ 4 1000 1000 4096 Oct 12 08:25 RS2_OK5491_PK71075_DK68880_FQ21_20090618_172447_HH_VV_HV_VH_SLC
drwx------ 4 1000 1000 4096 Oct 12 08:29 RS2_OK5491_PK71077_DK68882_FQ21_20090805_172450_HH_VV_HV_VH_SLC
drwx------ 4 1000 1000 4096 Oct 12 08:30 RS2_OK5491_PK71078_DK68883_FQ21_20090829_172451_HH_VV_HV_VH_SLC
drwx------ 4 1000 1000 4096 Oct 12 08:32 RS2_OK5491_PK71079_DK68884_FQ21_20091016_172455_HH_VV_HV_VH_SLC
drwx------ 4 1000 1000 4096 Oct 12 08:33 RS2_OK5491_PK86811_DK84561_FQ21_20090712_172447_HH_VV_HV_VH_SLC

The images are level one SLC (single look complex) format. For example, below are listed the contents of the image directory corresponding to acquistion date 20100426 (April 26, 2010). The four polarization combinations HH, HV,VH and VV are are stored as complex numbers in GeoTiff format:

In [3]:
ls -l /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC/
total 336276
-rw------- 1 1000 1000   741610 Apr 27  2010 BrowseImage.tif
-rw------- 1 1000 1000    44218 Oct  2  2009 DFAIT_RS2 EULA_Single User License.pdf
-rw------- 1 1000 1000      623 Jun 15 16:12 GEARTH_POLY.kml
-rw------- 1 1000 1000    89600 Apr 28  2010 Thumbs.db
-rw------- 1 1000 1000 85731468 Apr 27  2010 imagery_HH.tif
-rw------- 1 1000 1000 85731468 Apr 27  2010 imagery_HV.tif
-rw------- 1 1000 1000 85731468 Apr 27  2010 imagery_VH.tif
-rw------- 1 1000 1000 85731468 Apr 27  2010 imagery_VV.tif
-rw------- 1 1000 1000    49475 Apr 27  2010 lutBeta.xml
-rw------- 1 1000 1000    49475 Apr 27  2010 lutGamma.xml
-rw------- 1 1000 1000    49475 Apr 27  2010 lutSigma.xml
drwx------ 3 1000 1000     4096 Oct 11 15:25 polsarpro/
-rw------- 1 1000 1000   119099 Apr 27  2010 product.xml
-rw------- 1 1000 1000   160314 Jun 15 16:12 product_header.txt
-rw------- 1 1000 1000    15152 Jun 15 16:12 product_lut.bin
-rw------- 1 1000 1000    49765 Jun 15 16:12 product_lut.txt
-rw------- 1 1000 1000      316 Mar 24  2010 readme.txt
drwx------ 2 1000 1000     4096 Sep 17  2014 schemas/

The subdirectory polsarpro/T3 contains the polarimetric coherency matrix elements generated from the polarization combinations. They were generated with the PolSARpro software provided as freeware by the European Space Agency (ESA). (See the discussion below in the Section on the processing chain.) The PolSARpro image files are in ENVI format:

In [13]:
ls -l /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC/polsarpro/T3
total 61476
-rw------- 1 1000 1000 6986432 Oct 11 19:19 T11.bin
-rw------- 1 1000 1000     245 Oct 12 08:02 T11.bin.hdr
-rw------- 1 1000 1000 6986432 Oct 11 19:19 T12_imag.bin
-rw------- 1 1000 1000     246 Oct 11 19:19 T12_imag.bin.hdr
-rw------- 1 1000 1000 6986432 Oct 11 19:19 T12_real.bin
-rw------- 1 1000 1000     241 Oct 11 19:19 T12_real.bin.hdr
-rw------- 1 1000 1000 6986432 Oct 11 19:19 T13_imag.bin
-rw------- 1 1000 1000     246 Oct 11 19:19 T13_imag.bin.hdr
-rw------- 1 1000 1000 6986432 Oct 11 19:19 T13_real.bin
-rw------- 1 1000 1000     241 Oct 11 19:19 T13_real.bin.hdr
-rw------- 1 1000 1000 6986432 Oct 11 19:19 T22.bin
-rw------- 1 1000 1000     241 Oct 11 19:19 T22.bin.hdr
-rw------- 1 1000 1000 6986432 Oct 11 19:19 T23_imag.bin
-rw------- 1 1000 1000     246 Oct 11 19:19 T23_imag.bin.hdr
-rw------- 1 1000 1000 6986432 Oct 11 19:19 T23_real.bin
-rw------- 1 1000 1000     241 Oct 11 19:19 T23_real.bin.hdr
-rw------- 1 1000 1000 6986432 Oct 11 19:19 T33.bin
-rw------- 1 1000 1000     241 Oct 11 19:19 T33.bin.hdr
-rw------- 1 1000 1000      96 Oct 11 19:19 config.txt
-rw------- 1 1000 1000   18773 Oct 11 19:19 metadata.xml

We'll return to these images later in the tutorial. But first a little theory:

Statistical Properties of polSAR Images

Vector and matrix representations

A fully polarimetric SAR measures a $2\times 2$ scattering matrix $S$ at each resolution cell on the ground. The scattering matrix relates the incident and the backscattered electric fields $E^i$ and $E^b$ according to

$$ \pmatrix{E_h^b \cr E_v^b} =\pmatrix{S_{hh} & S_{hv}\cr S_{vh} & S_{vv}}\pmatrix{E_h^i \cr E_v^i}. $$

Here $E_h^{i(b)}$ and $E_v^{i(b)}$ denote the horizontal and vertical components of the incident (backscattered) oscillating electric fields directly at the target. These can be deduced from the transmitted and received radar signals via the so-called far field approximations. If both horizontally and vertically polarized radar pulses are emitted and discriminated then they determine, from the above Equation, the four complex scattering matrix elements. The per-pixel polarimetric information in the scattering matrix $S$, under the assumption of reciprocity ($S_{hv} = S_{vh}$), can then be expressed as a three-component complex vector

$$ s = \pmatrix{S_{hh}\cr \sqrt{2}S_{hv}\cr S_{vv}}, $$

where the $\sqrt{2}$ ensures that the total intensity (received signal power) is consistent. It is essentially these vectors which are provided in the SLC level one data discussed above. The total intensity is referred to as the span and is the complex inner product of the vector $s$,

$$ {\rm span} = s^\top s = |S_{hh}|^2 + 2|S_{hv}|^2 + |S_{vv}|^2. $$

This is a real number and the corresponding gray-scale image is called the span image. The observation vector $s$ can be shown to be a realization of a complex multivariate normal random variable. An equivalent and often preferred representation is in terms of the coherency vector

$$ k = {1\over\sqrt{2}}\pmatrix{S_{hh} + S_{vv}\cr S_{hh} - S_{vv} \cr 2S_{hv}}. $$

The polarimetric signal is can also be represented by taking the complex outer product of $s$ with itself:

$$ C = s s^\top = \pmatrix{ |S_{hh}|^2 & \sqrt{2}S_{hh}S_{hv}^* & S_{hh}S_{vv}^* \cr \sqrt{2}S_{hv}S_{hh}^* & 2|S_{hv}|^2 & \sqrt{2}S_{hv}S_{vv}^* \cr S_{vv}S_{hh}^* & \sqrt{2}S_{vv}S_{hv}^* & |S_{vv}|^2 }. $$

The diagonal elements of $C$ are real numbers, with span $= {\rm tr}(C)$, and the off-diagonal elements are complex. This matrix representation contains all of the information in the polarized signal.

Multi-looking

The matrix $C$ can be averaged over the number of looks (number of adjacent cells used to average out the effect of speckle) to give an estimate of the covariance matrix of each multi-look pixel:

$$ \bar{C} ={1\over m}\sum_{\nu=1}^m s(\nu) s(\nu)^\top = \langle s s^\top \rangle = \pmatrix{ \langle |S_{hh}|^2\rangle & \langle\sqrt{2}S_{hh}S_{hv}^*\rangle & \langle S_{hh}S_{vv}^*\rangle \cr \langle\sqrt{2} S_{hv}S_{hh}^*\rangle & \langle 2|S_{hv}|^2\rangle & \langle\sqrt{2}S_{hv}S_{vv}^*\rangle \cr \langle S_{vv}S_{hh}^*\rangle & \langle\sqrt{2}S_{vv}S_{hv}^*\rangle & \langle |S_{vv}|^2\rangle }, $$

where $m$ is the number of looks. This matrix (or rather the equivalent coherency matrix $\langle k k^\top \rangle$) is generated by PolSARpro and stored in the subdirectory polsarpto/T3 which we listed in a previous cell above. Rewriting the first of the above equalities,

$$ m\bar{C} = \sum_{\nu=1}^m s(\nu) s(\nu)^\top =: x. $$

This quantity $x$ is a realization of a complex sample matrix which is known to have a complex Wishart distribution with $N\times N$ covariance matrix $\Sigma$ (here $N=3$) and $m$ degrees of freedom:

$$ p_{W_c}( x) ={|x|^{(m-N)}\exp(-{\rm tr}(\Sigma^{-1} x)) \over \pi^{N(N-1)/2}|\Sigma|^{m}\prod_{i=1}^N\Gamma(m+1-i)},\quad m \ge N, $$

provided that the vectors $s(\nu)$ are independent and drawn from a complex multivariate normal distribution.

Equivalent number of looks

However when multi-look averaging takes place, the observation vectors $s(\nu)$ are not completely independent and will generally be correlated somewhat. In order to account for this, the complex Wishart distribution is often parameterized with ENL (rather than $m$) degrees of freedom, where ENL is the so-called equivalent number of looks. This quantity can be estimated from the image itself.

Dual and single polarimetric imagery

The scattering vector given above corresponds to so-called full, or quad polarimetric SAR. Satellite-based SAR sensors often operate in reduced, power-saving polarization modes, emitting only one polarization and receiving two (dual polarization) or one (single polarization). The look-averaged covariance matrices are reduced in dimension correspondingly. For instance for dual polarization with horizontal transmission and horizontal and vertical reception,

$$ \bar{C} = \pmatrix{ \langle |S_{hh}|^2\rangle & \langle S_{hh}S_{hv}^*\rangle \cr \langle S_{hv}S_{hh}^*\rangle & \langle |S_{hv}|^2\rangle }, $$

and, for single polarization and horizontal transmission/reception, we get simply the intensity image

$$ \bar{I} = \langle |S_{hh}|^2\rangle \quad {\rm or} \quad \bar{I} = \langle |S_{vv}|^2\rangle. $$

In the the dual pol case, the observations are complex Wishart distributed with $N=2$, in the single pol case they are gamma distributed.

Change Detection

Bitemporal data

The following is discussion is based on Conradsen et al (2004).

As we have seen, we can represent a pixel vector in an $m$ look-averaged polSAR image in covariance matrix format by $\bar C$, where

$$ m\bar C = x = \sum_{\nu=1}^m s(\nu) s(\nu)^\top $$

is a realization of a random matrix $X$ with a complex Wishart distribution.

In order to derive a change detection procedure for two polarimetric SAR images, we formulate a statistical test. For each pixel in the two image, define the null (or no-change) simple hypothesis as

$$ H_0:\quad \Sigma_1 = \Sigma_2 = \Sigma, $$

and the alternative composite hypothesis as

$$ H_1:\quad \Sigma_1 \ne \Sigma_2. $$

Under $H_0$ the maximum likelihood for $\Sigma$ can be shown to be given by

$$ L(\hat\Sigma) = { |x_1|^{m-3}|x_2|^{m-3}\exp(-2m\cdot{\rm tr}(I)) \over \left({1\over 2m}\right)^{3\cdot 2m}| x_1+ x_2|^{2m}\Gamma_3(m)^2 }, $$

where $I$ is the $3\times 3$ identity matrix and ${\rm tr}(I)=3$. Under $H_1$ the maximum likelihood for $\Sigma_1$ and $\Sigma_2$ is

$$ L(\hat\Sigma_1,\hat\Sigma_2) = { |x_1|^{m-3}|x_2|^{m-3}\exp(-2m\cdot{\rm tr}(I)) \over \left({1\over m}\right)^{3m}\left({1\over m}\right)^{3m} |x_1|^m |x_2|^m\Gamma_3(m)^2 } $$

Then the likelihood ratio test has the critical region for rejection of the no-change hypothesis

$$ Q = {L(\hat\Sigma) \over L(\hat\Sigma_1,\hat\Sigma_2) } = 2^{6m}{ |x_1|^m |x_2|^m \over |x_1 + x_2|^{2m} } \le k. $$

Finally, one can derive (Conradsen et al (2004)) the following approximation for the statistical distribution of the test statistic $Q$:

$$ {\rm prob}(-2\rho\log Q\le z) \simeq P_{\chi^2;N^2}(z) + \omega_2\left[ P_{\chi^2;N^2+4}(z) - P_{\chi^2;N^2}(z) \right], $$

where $P_{\chi^2;m}(z)$ is the chi square distribution wth m degrees of freedom, $$ \rho = 1 - {2N^2-1\over 6N}\cdot{3\over 2 m}, $$

and

$$ \omega_2 = - {N^2\over 4}\cdot\left(1-{1\over \rho}\right)^2 + {N^2(N^2-1)\over 24 \rho^2}\cdot{7\over 4m^2}. $$

In practice we choose a significance level $\alpha$, e.g., $\alpha = 0.01$, and decision threshold $z$ such that

$$ {\rm prob}(-2\rho\log Q\le z) = 1-\alpha $$

and interpret all pixels with larger values of $-2\rho\log Q$ as change.

Multitemporal data

The preceding discussion generalizes in a straightforward way to a time series of $k$ images (Conradsen et al, Determining the points of change in time series of polarimetric SAR data, to be published). In the equation for the test statistic $Q$ above, the numerator consists of a product $k$ determinants $|x_1|^m|x_2|^m\dots|x_k|^m$ and the denominator similarly the determinant of the sum of the $k$ observations $|x_1+x_2+\dots x_k|^{2m}$. The multitemporal test is referred to as the omnibus test will in general be more powerful (have a higher detection probability) for the same significance level than a bitemporal test just involving the first and last images in the series.

A Processing Chain

The change detection method implies the following processing sequence in order to generate a change map from a time series of polarimetric SAR images provided at the single look complex (SLC) processing level:

  1. First of all the multi-look polarimetric SAR images in covariance or coherency matrix format are generated from the SLC data with PolSARpro (or, if prefered, from the Sentinel 1 Toolbox also available from the European Space agency). Presently this must be done outside of the Docker container (and IPython) since the ESA software is only available in the form of a graphics interface. The coherenecy matrix has the same eigenvalues and hence the same determinant as the covariance matrix, so that the hypothesis test described above can be applied unchanged to either format. The rest of the processing takes place in the IPython notebook.

  2. The matrix images are imported by MapReady for georeferencing. The bash script /sar/mapready.sh in the main automates the procedure. MapReady will output the geocoded covariance/coherency matrix image in the form of co-registered GeoTiff files, one for each diagonal matrix element and two (real and imaginary parts) for each off-diagonal component. A python script /sar/ingest.py is called automatically to combine these files to a single multi-band image in floating point format.

  3. The ENL (equivalent number of looks) can (optionally) be estimated with the script enlml.py. A multivariate estimator is used as described by Anfinsen et al. (2009).

  4. The bitemporal change detection algorithm is invoked with the bash script /sar/wishart.sh. This script calls the Python programs /sar/register.py to co-register the two images and then /sar/wishart.py to perform the pixel-wise hypothesis tests. The test statistics $-2\rho\log Q$ and change probabilities ${\rm prob}(-2\rho\log Q\le z)$ are written to a two-band GeoTiff file. Additionally a change map showing changes at significance level 0.01 in red overlayed onto the span image is writen to a three-band (RGB) GeoTiff file. The multitemporal algorithm is invoked similarly with /sar/omnibus.sh, which calls /sar/register.py to co-register all of the images with the first in the series, and then /sar/omnibus.py to perform the multitemporal mnibus test.

An Example: Radarsat-2 Quadpol Images

Returning now to the Radarsat-2 image acquired April 24, 2010, we will geocode it with MapReady (step 2 in the above processing chain). The bash script /sar/mapready takes two arguments, the acquisition date in the format yyymmdd and the sensor (presently rs2quad or tsxdual):

In [1]:
!./mapready.sh 20100426 rs2quad
Geocoding polSARpro multilook polarimetric matrix image with mapready ...
Original SLC image dimensions:      rows 5539  cols 3788
After multi-looking with polSARpro: rows 1384  cols 1262
Azimuth looks: 4
Range looks:   3
***** processing polSARpro polarimetric matrix image:
***** /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC
***** ...
***** Done, see mapready.log
***** Combining into a single image file ...
=========================
     Ingest SAR
=========================
Wed Jun 17 08:50:17 2015
Directory /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/
writing band 1
writing band 2
writing band 3
writing band 4
writing band 5
writing band 6
writing band 7
writing band 8
writing band 9
elapsed time: 80.1424710751
Multiband image is /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif

We see that the multi-look images were created from PolSARpro data with $4\times 3 = 12$ looks. This corresponds to a square pixel size of close to $20.5\times 20.5$ meters. The combined coherency matrix image at this resolution is stored in geotiff format.

Before we can display the geocoded image, we have to enable Matplotlib functionality within the notebook with the so-called magic command

In [14]:
%matplotlib inline

We will use the Python script /sar/dispms.py for displaying. Here is the help:

In [15]:
run dispms -h
Usage: python dispms.py [-c] [-C] [-f filename1] [-F filename2] [-p posf] [P posF [-d dimsf] [-D dimsF]

                                        [-e enhancementf] [-E enhancement]

            if -f is not specified it will be queried

            use -c or -C for classification image
 
            RGB bandPositions and spatialDimensions are lists, e.g., -p [1,4,3] -d [0,0,400,400] 

            enhancements: 1=linear255 2=linear 3=linear2pc 4=equalization 5=logarithmic

The command below will generate an RGB color composite of the three diagonal elements:

In [4]:
run /sar/dispms -p [1,6,9] -f /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3//polSAR.tif

The scene above was acquired over the city of Bonn, Germany (upper right hand corner with the Rhine river). The gray, featurless areas are mixed forest.

To check the number of looks, We will run the /sar/enlml.py script (step 3 in the processing chain) on a spatial subset which includes a lot forested land cover. The spatial subset is entered with the -d flag as in the /sar/dispms script.

Here we choose -d [800,400,500,500]:

In [5]:
run enlml -d [800,400,500,500] /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3//polSAR.tif
=========================
     ENL Estimation
=========================
Tue Jun 16 07:44:11 2015
infile:  /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3//polSAR.tif
Quad polarimetry
filtering...
row:  100  200  300  400 
Mode: 5.6913998909
ENL image written to: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR_enl.tif
elapsed time: 71.0568621159

There are two modes (maxima) at about 6 and 12 looks. Here is the ENL image

In [8]:
run dispms -e 2 -f /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR_enl.tif

The bright areas correspond to mixed forest (around 12 equivalent looks) and allow for the best estimate of the ENL, see Anfinsen et al. (2009). So we can conclude that the nominal and eqivalent number of looks are the same and equal to 12.

Next we geocode two more Radarsat-2 images, the ones acquired on May 20, 2010 (20100520) and June 7, 2010 (20100707):

In [3]:
!./mapready.sh 20100520 rs2quad
Geocoding polSARpro multilook polarimetric matrix image with mapready ...
Original SLC image dimensions:      rows 5538  cols 3788
After multi-looking with polSARpro: rows 1384  cols 1262
Azimuth looks: 4
Range looks:   3
***** processing polSARpro polarimetric matrix image:
***** /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100520_172458_HH_VV_HV_VH_SLC
***** ...
***** Done, see mapready.log
***** Combining into a single image file ...
=========================
     Ingest SAR
=========================
Fri Jun 19 09:31:55 2015
Directory /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100520_172458_HH_VV_HV_VH_SLC_MapReady/T3/
writing band 1
writing band 2
writing band 3
writing band 4
writing band 5
writing band 6
writing band 7
writing band 8
writing band 9
elapsed time: 31.0119628906
Multiband image is /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100520_172458_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif
In [16]:
!./mapready.sh 20100707 rs2quad
Geocoding polSARpro multilook polarimetric matrix image with mapready ...
dos2unix: converting file /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100707_172459_HH_VV_HV_VH_SLC/polsarpro/T3/config.txt to Unix format ...
Original SLC image dimensions:      rows 5537  cols 3788
After multi-looking with polSARpro: rows 1384  cols 1262
Azimuth looks: 4
Range looks:   3
***** processing polSARpro polarimetric matrix image:
***** /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100707_172459_HH_VV_HV_VH_SLC
***** ...
***** Done, see mapready.log
***** Combining into a single image file ...
Overwriting previous combined image
=========================
     Ingest SAR
=========================
Thu Oct 15 14:18:33 2015
Directory /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100707_172459_HH_VV_HV_VH_SLC_MapReady/T3/
writing band 1
writing band 2
writing band 3
writing band 4
writing band 5
writing band 6
writing band 7
writing band 8
writing band 9
elapsed time: 1.36827516556
Multiband image is /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100707_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif

Now let's perform the last processing step, polSAR change detection. The bitemporal bash script /sar/wishart.sh needs five input parameters, the two acquistion times in yyymmdd format, a spatial subset, and the two ENL values (which in principle can be different). We choose first of all the April and June images:

In [17]:
!./wishart.sh 20100426 20100707 [400,400,1000,1000] 12.0 12.0
***** Bitemporal PolSAR Change Detection
 
***** registering ...
=========================
     Register SAR
=========================
Thu Oct 15 14:32:30 2015
Reference image: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif
Target image: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100707_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif
warping 9 bands (quad pol)...
elapsed time: 14.1322660446
Warped image written to: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100707_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR_warp.tif
***** complex Wishart change detection ...
============================================
Bi-temporal Complex Wishart Change Detection
============================================
Thu Oct 15 14:32:45 2015
first filename:  /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif
number of looks: 12.000000
second filename:  /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100707_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR_warp.tif
number of looks: 12.000000
Quad polarimetry
test statistic and change probabilities written to: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/wishart(20100426-20100707).tif
change map image written to: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/wishart(20100426-20100707)_cmap.tif
elapsed time: 2.00982999802

Here is the 1% significance level change map image generated by the above script:

In [19]:
run dispms -e 4 -p [1,2,3] -f /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/wishart(20100426-20100707)_cmap.tif

Over the interval late April - early June separating the two acquisitions the significant changes are mostly in the agricultural areas near center right and upper left. Barge movements on the Rhine river are clearly evident. Zooming in on the upper left hand corner we can see a flooded sand quarry pit with two dredging arms that are in continual motion, giving rise to significant change signals.

In [21]:
run dispms -e 4 -p [1,2,3] -d [1,1,400,400] -f /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/wishart(20100426-20100707)_cmap.tif

Finally, we'll run the omnibus test statistic on all three images. The script requires the acquisition times, the spatial subset, the number of looks and the desired significance level:

In [23]:
!./omnibus.sh 20100426 20100520 20100707 [400,400,1000,1000] 12 0.01
***** Multitemporal PolSAR Change Detection **********
******************************************************
number of images  3
ENL               12
spatial subset    [400,400,1000,1000]
=========================
     Register SAR
=========================
Thu Oct 15 14:42:29 2015
Reference image: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif
Target image: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100520_172458_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif
warping 9 bands (quad pol)...
elapsed time: 14.3113701344
Warped image written to: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100520_172458_HH_VV_HV_VH_SLC_MapReady/T3/polSAR_warp.tif
=========================
     Register SAR
=========================
Thu Oct 15 14:42:43 2015
Reference image: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif
Target image: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100707_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif
warping 9 bands (quad pol)...
elapsed time: 14.4490599632
Warped image written to: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100707_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR_warp.tif
===============================================
Multi-temporal Complex Wishart Change Detection
===============================================
Thu Oct 15 14:42:58 2015
first (reference) filename:  /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/polSAR.tif
number of looks: 12.000000
test statistic and change probabilities written to: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/omnibus(20100426-1-20100707).tif
elapsed time: 2.0478079319
change map image written to: /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/omnibus(20100426-1-20100707)_cmap.tif

We will display the bitemporal and multitemporal change maps side-by-side:

In [25]:
run dispms -e 4 -p [1,2,3] -E 4 -P [1,2,3] -d [1,1,400,400] -D [1,1,400,400] -f /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/wishart(20100426-20100707)_cmap.tif -F /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/omnibus(20100426-1-20100707)_cmap.tif 

The multitemporal test (right hand image) detects more changes than the bitemporal test, at the same 1% significance level.

In [ ]: